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ABSTRACT 


A  multilength  scale  method  based  on  asymptotic  expansion  homogenization  (AEH)  is  de¬ 
veloped  to  compute  minimum  energy  configurations  of  ensembles  of  atoms  at  the  fine  length 
scale  and  the  corresponding  mechanical  response  of  the  material  at  the  coarse  length  scale.  This 
multiscale  theory  explicitly  captures  heterogeneity  in  microscopic  atomic  motion  in  crystalline 
materials,  attributed,  for  example,  to  the  presence  of  various  point  and  line  lattice  defects.  The 
formulation  accounts  for  large  deformations  of  nominally  hyperelastic,  monocrystalline  solids. 
Unit  cell  calculations  are  performed  to  determine  minimum  energy  configurations  of  ensem¬ 
bles  of  atoms  of  body-centered  cubic  tungsten  in  the  presence  of  periodic  arrays  of  vacancies 
and  screw  dislocations  of  line  orientations  [111]  or  [100],  Results  of  the  theory  and  numerical 
implementation  are  verified  versus  molecular  statics  calculations  based  on  conjugate  gradi¬ 
ent  minimization  (CGM)  and  are  also  compared  with  predictions  from  the  local  Cauchy-Born 
rule.  For  vacancy  defects,  the  AEH  method  predicts  the  lowest  system  energy  among  the  three 
methods,  while  computed  energies  are  comparable  between  AEH  and  CGM  for  screw  dislo¬ 
cations.  Computed  strain  energies  and  defect  energies  (e.g.,  energies  arising  from  local  inter¬ 
nal  stresses  and  strains  near  defects)  are  used  to  construct  and  evaluate  continuum  energy 
functions  for  defective  crystals  parameterized  via  the  vacancy  density,  the  dislocation  density 
tensor,  and  the  generally  incompatible  lattice  deformation  gradient.  For  crystals  with  vacan¬ 
cies,  a  defect  energy  increasing  linearly  with  vacancy  density  and  applied  elastic  deformation 
is  suggested,  while  for  crystals  until  screw  dislocations,  a  defect  energy  linearly  dependent  on 
the  dislocation  density  tensor  appears  more  appropriate  than  the  quadratic  dependency  often 
encountered  in  the  continuum  plasticity  literature. 
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1.  INTRODUCTION 

Multiscale  methods  often  establish  constraints  that 
permit  numerical  juxtaposition  of  discrete  and  con¬ 
tinuum  material  descriptions.  This  necessarily  re¬ 
quires  enforcement  of  kinematic  approximations  to 
enforce  compatibility  between  two  otherwise  dis¬ 
parate  domains  to  reduce  the  overall  number  of 
degrees  of  freedom  and  produce  a  computation¬ 
ally  tractable  problem.  One  such  kinematic  approx¬ 
imation  is  the  Cauchy-Born  rule  (CB),  which  en¬ 
forces  tangent  maps  at  the  fine  scale  to  simulate 
the  locally  homogeneous  deformation  of  bulk  three- 
dimensional  crystals  [1,2].  Detailed  expositions  of 
the  underlying  theory  can  be  found  in  several  ref¬ 
erences  [3-5].  The  utility  of  the  CB  is  not  confined 
exclusively  to  bulk  crystals.  Following  extension  of 
the  Cartesian  theory  commonly  used  for  bulk  lat¬ 
tices  to  curvilinear  manifolds,  deformations  of  pla¬ 
nar  and  lower  dimensional  monocrystals  can  be  ef¬ 
ficiently  described  as  well  [6-8]. 

The  importance  of  these  types  of  kinematic  ap¬ 
proximations  manifests  in  their  continued  founda¬ 
tional  presence  in  multiscale  methods  for  materi¬ 
als  modeling.  Early  models  used  to  study  atom¬ 
istic  processes  in  the  vicinity  of  crack  tips  or  disloca¬ 
tion  cores  relied  heavily  on  continuum  elasticity  the¬ 
ory,  without  including  atoms  in  the  far-field.  Early 
work  featured  one-way  [9-14]  or  two-way  [15-19] 
coupled  methods,  in  which  displacement  fields  es¬ 
tablished  at  the  interface  between  continuum  and 
atomistic  regions  were  computed  either  from  so¬ 
phisticated  interfacial  conditions  or  from  initial  con¬ 
ditions  derived  from  continuum  elasticity  theory. 
Increases  in  computing  power  permitted  more  re¬ 
alistic  two-way  couplings,  whereby  atomistic  fields 
were  permitted  to  affect  the  far-field  elastic  continua 
through  the  latter's  discretization  with  finite  ele¬ 
ments  [20-23].  Such  improvements  in  the  coupling 
algorithms  enabled  description  of  dynamic  crack 
growth  [21].  These  approximations  were  also  pru¬ 
dent  at  the  time  as  uninteresting  effects  in  nearly 
homogeneously  deforming  regions  far  from  defects 
could  be  disregarded  and  because  of  the  limitations 
of  then  available  computing  resources. 

Novel  methods  involving  electronic  structure 
have  permitted  consideration  of  higher  accuracy 
calculations  in  the  near  field  [24,25],  with  the  finite 
element  domain  in  the  far-field  remaining  essen¬ 
tially  linear  elastic.  Although  some  early  work  at¬ 


tempted  to  parameterize  stress-strain  relationships 
in  the  far-field  region  with  atomic  potentials  [20- 
22],  more  efficient  considerations  of  selected  atom¬ 
istic  effects  on  material  behavior  were  achieved 
through  the  initial  developments  of  the  quasicon¬ 
tinuum  theory  [26-28]  that  employed  hyperelastic 
constitutive  behavior,  derived  from  atomistic  poten¬ 
tials,  for  the  overlaying  finite  elements.  In  the  sub¬ 
sequent  decade,  significant  new  developments  in 
methodologies  have  improved  the  fidelity  of  atom- 
istically  informed,  continuum  multiscale  computa¬ 
tional  methods  [29-31]. 

The  multiscale  modeling  approaches  discussed 
thus  far  reduce  the  problem  dimensionality  through 
various  kinematic  approximations  such  as  the  CB. 
Even  with  presumable  improvements  in  the  perfor¬ 
mance  of  future  generations  of  microprocessors,  the 
increasing  demand  of  higher-fidelity  modeling  [32] 
will  require  increases  in  computational  resources 
that  will  outpace  projected  computer  hardware  im¬ 
provements.  Such  demands  in  fidelity  mandate  im¬ 
provements  over  the  CB  in  situations  involving,  for 
example,  large  defect  densities  [33]  and  complex  ac¬ 
tive  lattices  [34]. 

Tewary  and  colleagues  [35,36]  demonstrated  the 
feasibility  of  representing  lattice  defects  in  reduced 
form  using  Green's  function  methods.  In  such  a 
lower-order  method,  the  defect  core  can  be  repre¬ 
sented  with  semianalytical  functions,  thereby  pro¬ 
viding  an  initial  estimate  that  can  be  used  as  an 
initial  state  in  a  more  computationally  intensive, 
larger-scale  atomistic  simulation.  However,  gen¬ 
eral  situations  involving  heterogeneous  strain  fields 
and  high  defect  densities,  as  opposed  to  isolated 
defects,  are  usually  not  considered.  Furthermore, 
increasing  evidence  supports  the  idea  that  concur¬ 
rent  simulations  at  finite  temperature  in  which  a  sta¬ 
tistical  (quantum  mechanical)  or  discrete  (molecu¬ 
lar)  domain  is  interfaced  with  a  continuum  domain 
presents  profound  obstacles  to  computer  method 
development  [37-40].  Although  the  present  investi¬ 
gation  is  restricted  to  quasi-static  isothermal  condi¬ 
tions,  it  is  noteworthy  that  self-consistent  hierarchi¬ 
cal  approaches  provide  a  means  of  preserving  sta¬ 
tistical  ensembles  of  atomic  motions  in  the  contin¬ 
uum  domain,  such  as  at  a  finite  element  quadrature 
point,  through  the  use  of  unit  cell-based  averaging. 

The  asymptotic  expansion  homogenization 
(AEH)  method  is  formulated  and  implemented 
in  the  present  article,  building  on  the  original 
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theory  of  Chung  and  colleagues  [41,42].  Discrete 
simulations  are  executed  at  the  atomistic  level,  with 
each  volume  element  (unit  cell)  of  atoms  subjected 
to  periodic  boundary  conditions.  Asymptotic  ho¬ 
mogenization  methods  [43,44]  are  used  to  compute 
the  macroscopic  tangent  stiffness  associated  with 
the  mechanical  response  of  the  ensemble  of  atoms 
in  the  unit  cell.  The  CB  kinematic  rule  is  invoked 
for  imposition  of  the  bulk  continuum  deforma¬ 
tion,  with  the  fine-scale  displacements  of  individual 
atoms  identified  with  the  inner  displacements  in  the 
asymptotic  approximation.  The  present  approach 
appears  ideal  in  its  present  form  for  addressing  the 
response  of  microstructures  containing  periodically 
distributed  defects,  in  contrast  to  other  methods 
(e.g.,  [26])  initially  developed  to  address  isolated 
defects.  This  is  because  only  one  or  a  few  defects 
need  be  simulated  explicitly  at  the  atomistic  level 
within  the  context  of  the  periodicity  assumption 
invoked  in  our  homogenization  scheme.  However, 
owing  to  this  very  same  periodicity  assumption, 
the  method  suffers  in  the  sense  that  isolated  (i.e., 
nonrepeating)  defects  cannot  be  easily  modeled. 
Specifically,  previous  applications  of  AEH  involved 
modeling  the  elastic  response  of  graphene  with 
point  defects  [41,42]  and  incorporation  of  the  kine¬ 
matics  and  energetics  of  finite  plastic  deformation 
in  metallic  crystals  at  the  coarse  scale  [45].  Other 
developments  have  appeared  in  the  generalized 
mathematical  homogenization  method  [46,47]  to 
more  carefully  account  for  the  dynamical  and 
thermal  behavior  of  atoms. 

Noteworthy  developments  discussed  in  the 
present  article  include  simulations,  in  deforming 
metallic  crystals,  of  vacancies  and  of  dislocations 
of  orientations  not  considered  previously  [45]  as 
well  as  parameterization  and  comparison  of  en¬ 
ergy  density  functions  from  continuum  defect  field 
(CDF)  theory  with  those  computed  using  the  AEH 
method.  Regarding  the  first  development,  the  AEH 
method  is  demonstrated  to  predict  atomic  configu¬ 
rations  of  unit  cells  containing  defects  with  nearly 
equal,  and  in  some  cases  lower,  system  energies 
than  conventional  lattice  statics  methods  [48,49] 
based  on  conjugate  gradient  minimization  (CGM). 
Comparison  with  results  from  CGM  is  of  interest 
because  CGM  is  often  invoked  at  the  fine  scale  in 
other  popular  multiscale  methods  such  as  the  quasi¬ 
continuum  theory  [26-28].  For  the  defects  consid¬ 
ered  here,  both  AEH  and  CGM  are  shown  to  predict 


more  realistic,  lower-energy  atomic  configurations 
than  the  local  CB  rule. 

Regarding  the  second  development,  some  back¬ 
ground  discussion  on  continuum  field  theories  of 
defects  is  now  warranted.  As  described  here,  these 
constitute  a  class  of  elastic-plastic  material  mod¬ 
els  for  crystalline  solids  that  are  also  multiscale 
theories,  though  discrete  atoms  are  not  involved 
at  the  fine  scale.  Rather,  effects  of  defects  and 
other  sources  of  microstructural  heterogeneity  are 
reflected  in  the  material  response  functions  (e.g., 
stored  energy,  yield  stress,  or  hardening  param¬ 
eters)  through  kinematic  variables,  internal  state 
variables,  and  spatial  gradients  of  these  variables 
[50-55].  Such  representations  account  for  strain 
hardening  in  plasticity  of  single  crystals  [52,53];  evo¬ 
lution  of  microstructures  associated  with  point,  line, 
and  surface  defects  [54,55];  and  regularization  of 
numerical  instabilities  in  continuum  inelasticity  im¬ 
plementations  [56,57], 

While  nonlinear  elastic  moduli  for  defect-free 
crystals  have  been  parameterized  versus  deforma¬ 
tion  [58,59],  continuum-level  energy  functions  of 
materials  with  periodically  distributed  defects  such 
as  vacancies  and  dislocations  have  heretofore  not 
been  explicitly  parameterized,  simultaneously,  ver¬ 
sus  defect  density  and  applied  deformation  using 
atomistic  modeling  tools.  In  many  cases,  sim¬ 
ple  phenomenological  formulations  are  used  for 
defect-dependent  strain  energy  in  continuum  mate¬ 
rial  models.  These  formulations  may  either  be  mo¬ 
tivated  from  microscopic  physics,  or  may  be  used 
simply  because  energetic  data  are  unavailable.  For 
porous  materials  (i.e.,  those  with  vacancy  defects), 
the  elastic  strain  energy  is  typically  reduced  linearly 
with  vacancy  concentration  such  that  stress  is  lin¬ 
early  reduced  [60]  along  the  lines  of  isotropic  dam¬ 
age  mechanics  theory  [61].  For  dislocations,  energy 
functions,  either  linear  or  quadratic,  in  the  line  den¬ 
sity  of  dislocations  per  unit  area  have  been  assumed. 
In  strain  gradient-based  continuum  dislocation  the¬ 
ories,  dislocation  populations  are  often  partitioned 
into  statistically  stored  dislocation  (SSD)  and  ge¬ 
ometrically  necessary  dislocation  (GND)  families, 
following  [62].  The  former  include  closed  dislo¬ 
cation  loops  and  contribute  no  net  Burgers  vector, 
while  the  latter  are  a  measure  of  the  incompatibil¬ 
ity  of  the  elastic  (or  plastic)  deformation  field.  Typi¬ 
cally,  a  linear  dependence  of  stored  (or  free)  energy 
of  the  crystalline  material  on  SSDs  or  total  disloca- 
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tion  density  [63]  is  assumed.  On  the  other  hand, 
quadratic  functions  are  often  used  to  account  for  the 
local  strain  energy  imparted  by  the  GND  tensor  [51- 
53,55].  This  form  is  used  for  simplicity,  motivated 
by  analogy  to  linear  elasticity  theory,  and  also  is 
used  to  provide  a  back  stress  dependent  on  the  den¬ 
sity  of  dislocations  [60,64]  or  its  gradient  [51,53].  Re¬ 
maining  open  issues  pertaining  to  continuum  mod¬ 
eling  of  dislocation  defects  in  the  context  of  gradi¬ 
ent  plasticity  are  quantification  of  the  length  scale 
parameter(s)  required  to  normalize  the  energy  [65] 
and  proper  selection  of  the  metric  tensor  used  to  col¬ 
lapse  the  GND  tensor  to  scalar  form  [55,66]. 

In  the  present  work,  numerical  results  obtained 
from  execution  of  the  AEH  multiscale  method  are 
used  to  motivate  the  choice  of  continuum  energy 
functions  and  associated  parameters,  in  particular 
for  single  crystalline  tungsten  (W)  containing  vacan¬ 
cies  or  screw  dislocations.  Only  energies  are  consid¬ 
ered,  and  not  stresses,  in  part  because  the  latter  are 
not  trivially  defined  from  an  atomistic  perspective 
[67].  The  investigation  is  limited  in  the  sense  that 
only  a  few  classes  of  defects  are  examined,  for  only 
one  material  (W),  and  these  defects  are  arranged  pe¬ 
riodically  in  the  lattice.  Furthermore,  all  simulations 
are  isothermal  at  null  temperature  (akin  to  molecu¬ 
lar  statics).  However,  this  work  constitutes  an  initial 
step  toward  computing  energies  used  in  continuum 
defect  theories  from  physics-based,  multiscale  com¬ 
putations,  as  opposed  to  phenomenological  curve 
fitting  of  the  material  response  to  macroscopic  stress 
data,  for  example.  Intermediate  scale  methods,  such 
as  phase  field  models  [68,69]  or  discrete  disloca¬ 
tion  simulations  [70-72],  may  ultimately  be  needed 
to  address  nonideal  defect  configurations  and  finite 
temperature  defect  kinetics  to  bridge  scales  of  atom¬ 
istic  resolution  and  continuum  crystal  mechanics  for 
arbitrarily  disordered  states  of  the  material. 

The  remainder  of  this  article  is  organized  as  fol¬ 
lows.  Section  2  features  derivations  of  the  multi¬ 
scale  homogenization  equations  applicable  to  geo¬ 
metrically  nonlinear  problems,  with  discrete  atoms 
resolved  at  the  fine  scale.  Section  3  includes  dis¬ 
cussion  of  numerical  implementation  and  special¬ 
ization  of  the  theory  to  tungsten  crystals.  Results 
of  demonstrative  simulations  of  defects  in  the  atom¬ 
istic  domain  are  given  in  Section  4.  These  results  are 
applied  toward  development  of  continuum-scale 
defect  field  descriptions  in  Section  5.  Vectors  and 
tensors  are  written  in  boldface  type,  with  scalars 


and  individual  components  of  vectors  and  tensors 
written  in  italic  font.  The  indicial  notation  is  fre¬ 
quently  employed,  with  summation  implied  over 
repeated  indices,  for  example,  AaBa  =  A 1 B\  + 
A2B2  +  A3Bs  when  a  =  1, 2, 3. 


2.  THEORY 


Reference  and  current  configurations  of  a  continu¬ 
ous  body,  denoted  by  Bq  and  B,  respectively,  are 
introduced.  Let  X  and  x  denote  coordinates  span¬ 
ning  the  reference  and  spatial  frames,  and  let  xa  = 
xa(XA,t)  denote  the  differentiable  motion  of  the 
material,  with  t  denoting  time.  The  deformation 
gradient  or  tangent  mapping  F  from  Bq  to  B  is  then 
written  as 


F  = 


<9x 

dX 


_  dx a_ 
■A  8XA 


(1) 


Strain  measures  are  introduced  as  follows,  where 
Gab  =  9aX  •  f)//X  and  gab  =  <9„x  •  dbx  are,  respec¬ 
tively,  metric  tensors  in  reference  and  spatial  coor¬ 
dinate  systems: 

Cab  =  FaAgabFbB  2EAB  =  FaAgabFbB  -  GAb  (2) 


Denoted  by  2,  P,  and  S,  the  Cauchy  stress,  first 
Piola-Kirchhoff  stress,  and  second  Piola-Kirchhoff 
stress,  respectively,  are  related  by 

=  j-lpo^pbA  =  j-lpaASABFbB  (3) 


Assuming  quasi-static  conditions,  local  forms  of  the 
balances  of  linear  and  angular  momentum  are  writ¬ 
ten  as  follows: 


P"A  +Ba  =  0  F"a  PbA  =  PaAFbA  (4) 

where  the  vertical  bar  denotes  covariant  differentia¬ 
tion  and  B  is  the  body  force  vector  per  unit  reference 
volume.  The  usual  symmetry  relations  for  stress 
tensors  hold  from  Eq.  (3)  and  the  second  of  Eq.  (4), 
£ab  =  £(ab)  and  SAB  =  S^AB\  where  parentheses 
indicate  symmetrization,  that  is,  ‘2A<ah>  =  Aab  +  Aba 
for  arbitrary  second-rank  tensor  A.  Multiplying  the 
first  of  Eq.  (4)  by  virtual  displacement  6u  and  inte¬ 
grating  over  reference  volume  V,  the  following  vir¬ 
tual  work  principle  is  obtained: 

J  PaBgab{8u)\BdV  =  J  Tagab8ubdA 

V  dV 

+  j  Ba gabbub dV  (5) 

V 
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with  the  traction  per  unit  reference  area  A  given  by 
Ta  =  PaB Nb,  where  NB  are  covariant  components 
of  the  unit  normal  vector  to  external  boundary  dV. 
Here  a  free  energy  potential  T  per  unit  reference 
volume  on  Bq  is  assumed  to  exist,  with  the  stress 
tensor  satisfying  the  following  hyperelastic  relation¬ 
ships: 


qAB  =  o  9*  =  ^ 

OCab  9Eab 


(6) 


For  the  particular  case  of  first-order  hyperelasticity, 
Eq.  (6)  becomes 


Next  an  additive  decomposition  of  displacements  at 
the  coarse  scale  is  assumed: 

/  =  r  +  5“  =  r  +  tva  (12) 

where  ua  represents  the  displacement  that  would 
exist  in  a  microscopically  homogeneous  medium 
and  ua  is  the  perturbation  in  displacement  due  to 
fine-scale  heterogeneity,  with  corresponding  fine- 
scale  representation  va.  The  corresponding  micro¬ 
scopic  decomposition  is 

va  =  £- v  =  va+va  =  (FaA  -  baA)  Ya  +  va  (13) 


SAB  = 


52T 


8Eab8Ecd 


Ecd  =  C 


ABCD 


Ecd  (7) 


where  CABCD  =  C (AB)(CD'>  is  the  fourth-rank  ten¬ 
sor  of  elastic  moduli  in  the  reference  frame.  Substi¬ 
tuting  Eq.  (6)  into  Eq.  (5)  and  using  Eq.  (3), 


/  2F°A  dCAB9ab  {bU)b\n  dV  =  J  Ta9ab^bdA 

V  dV 

+  j  Ba  gabbubdV  (8) 
v 


with  va  the  microscopic  displacement  arising  from 
the  projection  to  the  fine  scale  of  the  macroscopic  de¬ 
formation  gradient  FaA.  Differentiating  u  of  Eq.  (12) 
with  respect  to  XA  gives 


d 

8XA 


(ua  +  zva ) 


dua  |  8va 
dXA  +  Wa 


(14) 


where  we  have  appealed  to  the  second  of  Eq.  (9). 
The  left  side  of  Eq.  (8)  can  be  written  as  follows  in 
Cartesian  coordinates: 


v  v 


dT  d(& ua) 
dF%  dXB 


dV  (15) 


The  link  between  atomistic  (fine  scale)  and  con¬ 
tinuum  (coarse  scale)  resolutions  is  established  here 
via  the  AEH  technique  [41,42,45].  Let  fine  and 
coarse  length  scales  be  spanned  by  coordinates  ya  = 
ya{YA,t )  and  xa  =  xa(XA,t),  respectively.  Al¬ 
though  multiple  time  scales  have  been  used  else¬ 
where  [73],  in  the  present  scheme,  both  scales  are 
parameterized  by  the  same  temporal  variable  t. 
Multiscale  coordinates  are  related  by 

xa  =  tya  XA  =  eYA  (9) 

where  e  is  a  small  scalar  that  remains  constant 
throughout  the  time  history  of  deformation.  Coarse- 
and  fine-scale  displacements  u  and  v,  respectively, 
are  introduced.  These  are  restricted  below  to  coin¬ 
cident  Cartesian  coordinate  systems  in  the  reference 
and  spatial  frames: 

ua  =  xa  ~  baA XA  va=ya-  6 aAYA  (10) 

with  the  Cartesian  shifter  baA  =  1  for  a  =  A  and 
baA  =  0  for  a  A.  Corresponding  deformation  gra¬ 
dients  then  follow  as 

f)ua  f)ua 

+  ^  faA  =  wx  +  6^  <n> 


and  the  total  displacement  variation  6 ua  can  be  ex¬ 
pressed,  from  Eq.  (12),  as 

6  ua  =  bua  +  zbva  (16) 

Substituting  Eqs.  (15)  and  (16)  into  Eq.  (8)  produces 


v 


5T  d 
dFaB  dXB 


(■ bua  +  zbva)dV 


=  J  Tagab(bub  +  zbvb)dA 

dV 

+  j  Bagab(bub  +  tbvb)dV  (17) 

v 

which  is  then  volume  averaged  over  microdomain 
Y  to  yield 


1  f  f  chi'  f  d(bua)  ^  d(bva) 


dVdY 


Y  J  J  8FaB  \  dXB  8YB 

Y  V 

=  ^J  J  Tagab(bub  +  zbvb)dAdY 

Y  dV 

+  J  J  Ba9ab(t>ub  +  zbvh)dVdY  (18) 


Y  V 


Volume  5,  Number  3&4,  2007 


208 

Begell  House  Inc.,  http://begellhouse.com  Downloaded  2008-1-29  from  IP  128.63.66.91  by  Dr.  John  Clayton  (jdclayt) 


CHUNG  AND  CLAYTON 


Equation  (18)  is  satisfied  in  the  asymptotic  limit  e 
0  only  if 


y  v 


d^f  <9(6 ua) 
dFaB  dXB 


dVdY  =  lTagab8ubdA 


dv 


(19) 


+  /  Ba gabbubdV  (V8ua) 


v 


Y  V 


d'b  d(8va) 
dFaB  dYB 


dVdY  =  0  (V5va) 


(20) 


Solutions  of  Eqs.  (19)  and  (20)  converge  to  the  exact 
solution  (i.e.,  minimum  admissible  system  energy) 
for  ua  and  v"  when  the  displacement  field  is  peri¬ 
odic  in  Y,  for  example,  when  va  |  v-=0  =  v"  |  Y_  L  for 
a  domain  of  size  L.  Note  that  ua  is  constant  over  Y 
and  thus  automatically  satisfies  this  periodicity  con¬ 
straint. 

Equations  (l)-(20)  have  addressed  a  purely  con¬ 
tinuum  description  at  fine  and  coarse  length  scales. 
Presented  next  are  kinematic  and  thermodynamic 
assumptions  needed  to  relate  atomistic  and  contin¬ 
uum  fields  and  incrementally  update  atomic  coor¬ 
dinates.  Assume  that  in  reference  configuration  Bo, 
the  representative  volume  or  unit  cell  for  homoge¬ 
nization  consists  of  atoms  arranged  in  a  lattice,  per¬ 
haps  imperfect  due  to  the  presence  of  defects.  Fur¬ 
thermore,  assume  that  in  deformed  configuration  B, 
the  same  mass  and  number  of  atoms  exist  in  this 
representative  volume.  The  position  vector  for  each 
atom  j  in  configuration  Bo  is  given  by  7>(j)  =  aY  E *, 
where  angled  brackets  are  reserved  for  atomic  labels 
which  span  1  to  N,  a'Y  are  integers,  and  E,  are  the 
lattice  vectors  in  the  reference  configuration.  Spatial 
positions  of  atoms  in  configuration  B  are  then 
found  as  follows,  in  Cartesian  coordinates: 


*&>  =  S“a^>+«&>  (21) 

with  q/J;  being  a  displacement  vector  between  ref¬ 
erential  and  current  states  for  atom  j.  Let  and 

denote  vectors  separating  atoms  j  and  k  in  re¬ 
spective  configurations  Bo  and  B,  that  is, 

R0'\fc)  =  Z(fc)  “  zb>  (22) 

r(j\k)  =  z(k)  -  z<j)  (23) 

The  CB  rule  states  that  the  atoms  in  the  lattice  match 
the  gross  deformation  gradient  F  %  of  the  contin¬ 
uum,  ea  =  F  (( E  ' ,  with  e"  the  lattice  vectors  in 


the  deformed  configuration.  For  certain  crystals, 
namely,  defect-free  centrosymmetric  crystals,  this 
has  been  shown  to  be  an  effective  means  of  repre¬ 
senting  crystal  distortions  [4].  However,  symmetry- 
eliminating  or  symmetry-reducing  features,  such  as 
crystal  defects,  invalidate  CB  and  require  an  atom¬ 
istic  minimization  computation,  such  as  that  advo¬ 
cated  here  in  what  follows,  to  determine  distortions 
of  local  crystal  structure  near  defects. 

The  point  of  departure  here,  therefore,  is  to  use 
homogenization  to  determine  a  correction  to  the  lo¬ 
cal  region  near  the  defect  by  approximating  the  new 
coordinates  of  the  displaced  lattice  points  with 


z{o) 


=  F, 


(jk)A*{k) 


\o) 


(24) 


with  summation  implied  over  repeated  atomic  in¬ 
dices.  The  first  term  on  the  right-hand  side  of 
Eq.  (24),  A  Zf, ,  accounts  for  the  uniform  projec¬ 

tion  over  each  periodic  cell  of  the  macroscopic  lat¬ 
tice  deformation  field  to  the  fine  scale  (i.e.,  the  CB 
rule),  and  is  the  discrete  atomistic  analog  of  the 
perturbation  in  displacement  due  to  the  microscopic 
heterogeneity  given  previously  in  Eq.  (13),  written 
here  for  atom  j.  Spatial  separation  vectors  in  Eq.  (23) 
then  become 


rU\k)  =FR(i\fe)  +v<fc>  -  ■ vy)  =  FR(j\fc)  +?<j\fc>  (25) 

with  i(j\k)  =  accounting  for  deviations 

from  the  CB  rule  in  a  local  sense. 

Henceforth  in  the  present  work,  we  assume  a  free 
energy  potential,  measured  either  per  unit  reference 
volume  or  per  atom,  of  the  simple  form  depending 
only  on  the  relative  positions  of  atomic  nuclei  (Born- 
Oppenheimer  lattice  statics): 

M'  =  T'(q0->,z<j>) 

=  4'  (r(i\2),r<i\3)> r(2\3)>  ...r<N_i\N))  (26) 

Interatomic  forces  f(j)a  and  the  Hessian  matrix  (i.e., 
atomic  stiffness)  H/jk)ab  are  given  by 

d  T  <92T 

f{j)a  =  dg^  HWab  =  dq^dqbk)  (27) 

These  quantities  arise  directly  from  a  Taylor  series 
expansion  of  Eq.  (27)  about  a  fixed  set  of  reference 
coordinates  |  comprising  a  perfect  lattice: 

4'(q<j>)=4'o+/<fc>a|0  <f(k)+\ H(jk)ab\0  -  (28) 
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where  the  subscript  0  denotes  quantities  evaluated 
at  q/j)  =  0.  The  expansion  of  the  strain  energy  po¬ 
tential  of  the  solid  about  the  undeformed  state  is 
also  given  by  the  following  continuum  approxima¬ 
tion: 


1 

A 


dY  V 
1 

~  Y 


ST 

dF“A 


bva  N AdV dA 


d  f  ST 
8YA  I  dF\ 


y  v 


bvadVdY 


(33) 


T(n)  =  T0+  — 

1  S2T 
+  2  dFaAdFbB 


?>aA){FbB 


where  the  reference  energy  To  is  typically  assigned 
a  value  of  zero  in  strictly  continuum  models.  Note 
that  a  more  specific  objective  form  of  Eq.  (29)  could 
be  written  in  terms  of  the  strain  measure  Eab  of 
Eq.  (6),  for  example,  as  opposed  to  the  deformation 
gradient  FaA.  While  only  the  former  is  invariant  un¬ 
der  rigid  body  motion,  Eq.  (29)  suffices,  for  illustra¬ 
tive  purposes,  in  the  present  context.  On  taking  the 
first  variations  of  Eqs.  (28)  and  (29),  atomic  forces 
and  stresses  vanish  in  a  perfect  lattice  in  the  unde¬ 
formed  state,  that  is. 


Localizing  the  volume  integral  on  the  right-hand 
side  of  Eq.  (33)  and  considering  all  admissible  varia¬ 
tions  6v,  the  microscopic  linear  momentum  balance 
becomes 


d  /  8T  \ 
W A  \dFaA) 


dPaA 

dYA 


0  (inY) 


(34) 


as  the  area  integral  vanishes  since  one  may  select 
5v  =  0  on  A  =  dY.  Equation  (34)  is  general  in  the 
sense  that  no  assumption  is  made  on  the  order  of  the 
incremental  elastic  response;  for  example,  nonlinear 
higher-order  elastic  constants  are  admitted.  For  al¬ 
gorithmic  purposes,  however,  it  is  advantageous  to 
assume  a  first-order  hyperelastic  response  along  the 
lines  of  Eq.  (31),  an  appropriate  assumption  for  most 
engineering  metals  undergoing  quasi-static  defor¬ 
mations: 


hi)a\ o  —  0 


<9T 
d  FaA 


(30) 


leaving  the  following  force-displacement  and  stress- 
deformation  relations  on  neglecting  higher  than 
second-order  terms  in  Eqs.  (28)  and  (29): 


f(j)a  —  H(jk)ab 


0  q(k) 


jd-A _ (pAB 

ra  ~  ^ab 


(FbB-bbB)  (31) 

6  “a 


where  the  mixed-configurational  elastic  moduli  are 
CAB  =  <92T /dFaAdFhB.  Unobstructed  minimization 
of  Eqs.  (28)  and  (29)  may  proceed  only  when  H^k^ab 
and  CAB  are  positive  semidefinite  functions  of  their 
arguments.  Using  Eq.  (30)  and  equating  reference 
energies  T0  in  Eqs.  (28)  and  (29),  to  second  order. 


H{jk)abqa{j)q\k)  =  CAB  (FaA  -  baA )  ( FbB  -  bbB )  (32) 

a  relation  that  is  generalized  here  to  hold  when 
H{jk)ab  and  CAB  are  not  evaluated  at  the  reference 
state  (i.e.,  secant  moduli)  and  when  the  lattice  is  not 
initially  free  of  defects.  Integrating  the  right-hand 
side  of  Eq.  (15)  by  parts  and  applying  the  divergence 
theorem  over  volume  Y  with  oriented  surface  ele¬ 
ment  N Ad  A, 


-^A^abB{FbB~bbB))=  0  (35) 

where  CAB  is  a  mixed-variant  effective  elastic  mod¬ 
ulus  tensor.  Next,  from  Eq.  (11)  and  the  displace¬ 
ment  gradient  expression  Eq.  (14),  Eq.  (35)  becomes 


dCAbB  (  dub\  _d_  / FAB  dib  \ 
dYA  \dXB )  dYA  \  ab  3YB ) 


(36) 


The  intention  of  the  present  derivation  is  expres¬ 
sion  of  Eq.  (36)  in  terms  of  atomistic  and  macro¬ 
scopic  displacements.  Invoking  the  chain  rule  gives 


d  dYA  d  A  d 

dZB  ~  dZB  ~  U)BgyA 
u)  u> 


(37) 


where  the  linear  operator  n  simply  refers  initial 
position  ZBj  of  discrete  atom  j  to  fine  scale  refer¬ 
ence  coordinate  YA.  On  inserting  Eq.  (37)  into  the 
left-hand  side  of  Eq.  (36), 


dtAB  (  dub\  dtAB  (  dub  \ 
dYA  \8XB )  8ZA}  \dXB  J 

<92T  /  dub  \  _  B  /  dub  \ 

~  dqa{j)dFbB  ydxF ^ ab  ydXB  J 


(38) 
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where  ^fp)ah  =  —d2'l’/dq^dFbA.  On  appealing  to 
Eq.  (32),  and  assuming  small  atomic  perturbations 
v?m ) ,  the  right-hand  side  of  Eq.  (36)  becomes 

9  ( ipAB  9vb  \  .  9  ( tpAB  di,b(k)  \ 

8YA  \  ab  8YB J  dZfk)  \ab  dZf*}  j 

dH  ~b  _  ~b  root 

~  Qaa  a  b  V{k)  ~  H(kj)abV(k)  (39) 

°q(k)0q(j) 

In  Eqs.  (38)  and  (39),  the  — >  notation  denotes  the 
transformation  steps  vb  — >  vbk,  and  6  g 
Using  Eqs.  (38)  and  (39),  the  fine-scale  equilibrium 
Eq.  (36)  finally  becomes 

dub 

Dti)abgj[A  =  H(kj)abv\k)  (40) 

Subsequently,  Eq.  (40)  is  solved  for  the  inner  dis¬ 
placements  vbk,  describing  the  locally  perturbed 
atomic  coordinates.  Note  that  by  hypothesis  and 
analogy  with  the  two-scale  continuum  homogeniza¬ 
tion  theory,  solution  of  Eq.  (40)  converges  to  the  ex¬ 
act  solution  (i.e.,  lowest-energy  configuration  of  ad¬ 
missible  sets  of  atomic  coordinates)  only  when  vbkj 
is  periodic  over  Y.  In  other  words,  periodic  bound¬ 
ary  conditions  on  atomic  displacements  [74]  must 
be  applied  at  the  fine  scale  for  Eq.  (40)  to  yield  a 
meaningful  solution. 

Now  reconsider  coarse-scale  Eq.  (19),  assuming 
that  v  is  known  from  solution  of  Eq.  (40).  The  left- 
hand  side  of  Eq.  (19)  can  be  written  as 


(M*)  ^dvdY=. yT“»6“‘" 


+ j B‘s^ubdV+  yJJ, (43) 

V  YV 


3.  IMPLEMENTATION 

The  focus  of  the  present  effort  is  application  of 
Eq.  (40)  to  update  the  configuration  of  N  atoms 
subjected  to  macroscopic  displacement  gradient 
dub/dXA,  applied  uniformly  over  a  single  coarse- 
scale  integration  point,  such  that  Eq.  (43)  need  not 
be  solved  explicitly.  The  atoms  comprise  a  peri¬ 
odic  unit  cell  of  volume  Y  and  may  be  arranged 
initially  to  encompass  defects,  such  as  vacancies  or 
dislocations,  at  t  =  0.  Perturbations  from  homo¬ 
geneous  deformation  (i.e.,  perturbations  from  the 
CB)  are  measured  by  the  atomic  quantity  vb,^.  Note 
that  when  the  CB  is  accurate,  for  example,  when  de¬ 
fects  are  absent,  the  macroscopic  stress  at  a  given 
applied  deformation  level  is  minimized  with  respect 
to  variations  in  atomic  degrees  of  freedom  such  that 
d  (d^/dFaA)  /dqh(:i)  =  0,  and  thus  Eq.  (40)  yields 
vb{])  =  0.  ' 

Since  the  material's  mechanical  response  is  non¬ 
linear  in  the  presence  of  defects,  an  iterative  scheme 
is  employed  for  application  of  the  AEH  computa¬ 
tional  method  to  deforming  crystals.  Let  S  be  the 
enumerated  set  of  all  conceivable  nonuniform  de¬ 
formations  of  the  lattice: 


1 

Y 


Y  V 

1 

'  Y 


9T  5(6 ua) 

dF%  dXB 


dVdY 


(41) 


'iBA 

■'ab 


(  dub 

l  <9X^ 


dvb  \  <9(6 ua) 

W A  J  dXB 


dVdY 


Y  V 


where  fine-scale  atomic  coordinates  affect  FA  and 
CAB  through  Eqs.  (14)  and  (28),  respectively, 
and  thus  affect  the  macroscopic  response.  From 
Eqs.  (37)-(40), 


hi) 


=  f; 


’(jk)A^(k)  +  (Vi  €  S) 


(44) 


Equilibrium  Eq.  (40)  then  becomes 


Da  — 
{, j)ab  QXa 


H{jk)ab h  ^(fc) 


(45) 


where  the  line  search  parameter  q*  is  determined  it¬ 
eratively  such  that  a  local  minimum-energy  config¬ 
uration  is  attained: 


C 


BA 

ab 


dvb 

g^A 


C 


BA 


<92T 


dZ& 


dqfydF^ 


4> 


—  _  Tln 

U(:)abV(j) 


(42) 


The  coarse-scale,  static  linear  momentum  balance, 
Eq.  (19),  then  becomes 


min'T  (z  (q1))  ~  min'T  (z  (q/q))  (46) 

The  physical  implication  of  Eq.  (46)  is  that  the  min¬ 
imum  energy  state  corresponding  to  the  vector 
of  Eq.  (21)  can  be  approximated  by  the  state  that 
corresponds  to  Eq.  (44),  where  the  information  in 
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the  initial  solution  vector  vh/ky  determined  earlier  in 
the  computation  from  Eq.  (40),  has  been  reused.  A 
scenario  in  which  such  a  scheme  may  offer  signif¬ 
icant  computational  benefit  is  mechanical  calcula¬ 
tions  of  lattice  deformations  for  generating  stress- 
strain  curves,  whereby  the  deformation  gradient  is 
monotonically  increased. 

The  numerical  algorithm  proceeds  as  follows. 
First,  Eq.  (40)  provides  the  trajectory  for  v\k)  in  3  x  N 
solution  space.  Then  atomic  coordinates  z,y  and 
corresponding  displacements  q/q  are  updated  it¬ 
eratively  via  Eq.  (44).  With  each  iteration,  the  en¬ 
ergy  T  (z  (q/j\))  is  computed  using  the  updated 
coordinates  Z(j)  =  Z/(/  +  q^\ .  A  bisection  algo¬ 
rithm  is  used  to  efficiently  determine  the  particu¬ 
lar  value  of  q*  that  gives  the  minimum  energy  con¬ 
dition  dT  (z  (q*))  /dr)*  =  0  for  each  increment  in 
macroscopic  deformation  u. 

The  noteworthy  feature  of  the  present  AEH  com¬ 
putational  method  is  its  straightforward  trajectory, 
determined  by  Eq.  (45),  to  the  solution  of  updated 
atomic  coordinates  and  the  corresponding  energy 
minimum.  This  is  in  contrast  to  traditional  molecu¬ 
lar  dynamics  simulations  that  rely  on  a  series  of  in¬ 
cremental  loading  and  equilibration  steps  or  CGM 
methods  that  rely  on  series  of  incremental  straining 
and  minimization  steps,  both  in  full  3  x  N  solution 
space.  The  methods  are  compared  qualitatively  in 
Fig.  1.  For  the  molecular  methods,  the  stepwise  so¬ 
lution  path  denotes  successive  incremental  loading 
and  relaxation.  For  the  AEH  method,  the  trajec¬ 
tory  of  the  solution  is  denoted  by  the  correspond¬ 
ing  dotted  line,  with  distance  along  this  path  mea¬ 


strain  - ► 

FIGURE  1.  Schematic  comparison  of  numerical  ap¬ 
proaches 


sured  in  practice  by  the  value  of  line  search  parame¬ 
ter  r)\  It  is  also  important  to  note  that  no  assumption 
of  linearity  of  material  behavior  is  made.  Indeed, 
the  determination  of  the  numerical  value  q*  through 
the  condition  (z  (q*))  /<9q*  =  0  leaves  the  full 
material  nonlinear  (or  nonquadratic)  features  of  the 
atomistic  energy  T  undisturbed. 

The  AEH  method  is  applied  here  in  a  study  of 
the  mechanical  behavior  of  pure  tungsten  (W),  a 
BCC  transition  metal  of  relatively  high  mass  density. 
Its  combination  of  high  density,  high  strength,  and 
high  melting  point  render  it  a  popular  material  for 
use  in  defense  applications  such  as  ordnance  [63,75]. 
The  potential  energy  function  used  here  for  describ¬ 
ing  atomistic  interactions  in  W  is  discussed  in  what 
follows. 

An  empirical  N-body  potential  specifically  devel¬ 
oped  for  transition  metals  [76]  is  used  to  compute 
the  free  energy  T  of  Eq.  (26),  in  particular  because 
of  its  adequacy  for  describing  energies  and/ or  mo¬ 
tion  of  dislocations  and  other  lattice  defects  in  W,  as 
reported  elsewhere  [77-80].  According  to  the  rep¬ 
resentation  of  Finnis  and  Sinclair  [76],  the  summed 
total  potential  energy  E  of  a  set  of  atoms  at  positions 
{z(j)}  for  j  =  1,2,...,  His  given  by 


E  =  En  +  Ep  (47) 

where  fiy  is  the  N-body  term  that  is  a  function  of  the 
superposition  of  the  local  electronic  charge  densities 
Pif),  the  latter  obtained  from  summation  of  atomic 
charge  densities  4>.  Ep  is  a  pair  potential  that  ac¬ 
counts  for  core-core  interactions.  Specifically,  the  N- 
body  term  is 


EN  =  -Aj2f(PU))  (48) 

O') 


where 


/  (po>)  =  VPU)  P(j)  =  H  ^  (rovo) 

k 

is  always  nonnegative  and  real, 

r0\fc>  =  lr0\fc>l  =  lzw  -z0>l 


(49) 


(50) 


f  (r  —  d)2,  r  <  d 
\  0,  r  > 


(51) 


and  A  is  an  empirical  constant.  The  parameter  d  de¬ 
notes  an  adjustable  cutoff  zone  for  superposition  of 
local  charge  densities,  chosen  here  to  lie  between  the 
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second  and  third  nearest  neighbors,  that  is,  a  <  d  < 
a\J 2,  with  a  the  lattice  parameter.  The  pair  potential 
is  constructed  as 

Ep  =  \  Y  ^  (r<i\fc>)  <52) 

where  is  of  the  following  polynomial  form: 

rlt  (r)  =  /  (r  -  c)2  (Co  +  cir  +  C2r2)  ’  r  ^  c  (53) 
w  [  0,  r  >  c 


dr(M 

dF°A 


=  R 


(j\k)°.a 


dqb 


U) 


dF°A 


=  zw 


(58) 


Derivatives  of  E  with  respect  to  atomic  displace¬ 
ments  are  listed  in  [76]  and  are  not  repeated 
here.  In  the  present  numerical  scheme,  Eqs.  (55)- 
(58)  are  evaluated  analytically  and  subsequently 
used  in  fine-scale  equilibrium  Eq.  (45). 


4.  NUMERICAL  RESULTS  AND  VALIDATION 


with  empirical  constants  Co,  c\,  and  ru.  The  cutoff 
parameter  c,  like  d  of  Eq.  (51),  is  also  assigned  a 
value  between  the  second  and  third  nearest  neigh¬ 
bor  distances.  Finnis  and  Sinclair  [76]  determined 
the  other  constants  via  calibration  to  experimentally 
determined,  macroscopic  elastic  properties  for  sin¬ 
gle  crystalline  W.  Tables  1  and  2  list  the  experimen¬ 
tal  and  fitted  parameters,  respectively. 

Notice  that  E  is  the  total  energy  of  N  atoms  in 
the  fine-scale  unit  cell.  The  Helmholtz  free  energy 
density  in  a  continuum  sense,  T,  is  related  to  E  as 

*  =  U-^d  =  ^  <54) 

where  U  is  the  continuum  internal  energy,  r\  is  the 
continuum  entropy  per  unit  volume,  and  0  is  the  ab¬ 
solute  thermodynamic  temperature  of  the  system, 
which  we  assume  is  zero  in  the  last  of  Eq.  (54)  for 
the  present  implementation  in  the  context  of  molec¬ 
ular  statics,  such  that  E  =  [31417.  Also  included  in 
Eq.  (54)  is  the  scalar  |3  =  0.5a3,  a  constant  denot¬ 
ing  the  volume  occupied  by  each  atom  in  a  perfect 
reference  lattice.  This  unit  conversion  factor  was  in¬ 
cluded  implicitly  in  earlier  work  [45]. 

Substitution  of  Eqs.  (47)  and  (54)  into  earlier  def¬ 
initions  then  yields 
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(jk)ab 


J-J/A  \  „  h  - 


dq^dql) 
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(55) 


(56) 


(57) 


To  first  order,  from  Eqs.  (21)-(25),  with  7 A  ~ 
paAb(jk),  the  following  geometric  relationships 
arise: 


The  AEH  multiscale  computational  method  is  ap¬ 
plied  here  to  address  the  nonlinear  elastic  response 
of  body-centered  cubic  (BCC)  tungsten  (W)  con¬ 
taining  periodically  distributed  vacancies  and  screw 
dislocations  of  two  orientations.  In  these  simula¬ 
tions,  unit  cells  are  deformed  in  uniaxial  stretch  to 
2.5%  elongation.  The  primary  solution  variable  of 
interest  is  the  strain  energy  density  of  the  material 
in  the  presence  of  defects  contained  within  the  unit 
cell.  The  present  approach  readily  enables  paramet¬ 
ric  variations  of  the  defect  density  via  the  prescrip¬ 
tion  of  the  number  of  atoms  in  the  fine-scale  repre¬ 
sentation  relative  to  the  total  number  of  defects  em¬ 
bedded  within  the  unit  cell. 

One  intention  of  the  present  investigation,  dis¬ 
cussed  in  more  detail  in  Section  5,  is  examination 
of  aspects  of  material  behavior  that  could  be  used 
subsequently  in  stand-alone  continuum  defect  theo¬ 
ries,  in  particular,  details  associated  with  stored  en¬ 
ergy  of  defect  fields  and  the  effects  of  applied  defor¬ 
mations  on  elasticity  and  stored  energy  for  various 
fixed  defect  concentrations.  Despite  the  presump¬ 
tion  of  zero  temperature,  one  may  still  draw  conclu¬ 
sions,  at  least  in  a  qualitative  sense,  regarding  stress 
and  energetics  associated  with  dislocation  behavior 
in  BCC  metals  in  the  context  of  lattice  statics  calcula¬ 
tions  [74,78].  A  standard  trend  in  the  literature  has 
been  development  of  such  scaling  methods  for  the 
purely  mechanical  problem  before  extending  to  the 
finite  temperature  regime  [81]. 

Simulations  of  the  deformation  of  unit  cells  con¬ 
taining  various  numbers  of  atoms,  configured  to 
represent  several  classes  of  crystal  defects,  are  con¬ 
ducted.  In  the  calculations,  the  energy  density  and 
tangent  stiffness  of  a  single  Lagrangian  finite  ele¬ 
ment  integration  point  are  determined  by  the  mi¬ 
croscopic  (i.e.,  atomistic)  response.  Initial  atomic 
coordinates  are  found  using  a  two-step  procedure: 
first  the  linear-elastic  solution  for  displacement  field 
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TABLE  1.  Experimental  lattice  quantities  for  W 


a 

Lattice  parameter  [A] 

3.1652 

Eu 

Cohesive  energy  [eV/atom] 

-8.90 

Cu 

Elastic  constant  [GPa] 

522.4 

Cl  2 

Elastic  constant  [GPa] 

204.4 

C44 

Rhombohedral  shear  modulus  [GPa] 

160.6 

b 

Tetragonal  shear  modulus  [GPa] 

159.0 

B 

Bulk  modulus  [GPa] 

310.4 

Pu 

Cauchy  pressure  = (C12  —  C44)  [GPa] 

21.9 

TABLE  2.  Constants  for  Finnis-Sinclair  [76]  EAM  potential  (W) 


d[  A] 

4.40024 

A[eV] 

1.896373 

c  [  A] 

3.25 

Co 

47.1346499 

Cl 

-33.7665655 

C2 

6.2541999 

of  the  defect  is  applied  to  the  atoms,  then  a  conju¬ 
gate  gradient  algorithm  [49]  is  invoked  to  transition 
the  atomic  positions  to  a  stable  local  minimum  en¬ 
ergy  state.  Subsequently,  the  response  to  applied 
deformation  is  computed  using  our  AEH  scheme 
according  to  the  numerical  procedure  described  in 
Section  3.  The  applied  (i.e.,  coarse  scale)  deforma¬ 
tion  gradient  field  (in  conjunction  with  fine-scale 
periodicity)  is  uniaxial  stretching  over  a  range  of 
1.000  <  f  j  i  <  1.025,  with  the  lateral  edges  fixed  (co¬ 
variant  Cartesian  notation  is  used  here  and  in  sub¬ 
sequent  figures  for  simplicity,  i.e.,  F\  —>  Fn  ).  We 
also  compare,  for  validation  purposes,  the  final  con¬ 
figurations  attained  using  our  procedure  with  those 
obtained  from  incremental  energy  minimization  us¬ 
ing  the  CGM  method.  In  the  latter  approach,  a  small 
increment  in  the  stretch  field  is  first  imposed  uni¬ 
formly  over  all  atoms,  and  then  a  conjugate  gradi¬ 
ent  program  [49]  is  used  to  update  the  atomic  co¬ 
ordinates  to  the  corresponding  local  minimum  en¬ 
ergy  state.  This  process  continues,  with  a  new  set 
of  conjugate  gradient  minimization  iterations  con¬ 
ducted  on  application  of  each  successive  stretch  in¬ 
crement,  until  the  final,  fully  deformed  configura¬ 
tion  is  reached. 

Two  orientations  of  atomistic  unit  cells  are  inves¬ 
tigated.  In  the  first,  shown  in  Fig.  2,  the  axis  of  ap¬ 
plied  stretch  is  oriented  along  the  [11  Indirection  in 
the  lattice.  Here  the  BCC  unit  cell  is  a  rectangle  of  di¬ 
mensions  Li  x  L2  x  L3  =  a\/3Ni  x  C1V6N2  x 


where  N\,  N2,  and  A3  are,  respectively,  the  num¬ 
ber  of  repeating  planes  stacked  in  the  [111]-,  [112]-, 
and  [110]-directions.  In  the  second  orientation  (not 
shown),  the  axis  of  applied  stretch  is  oriented  along 
the  [100]-direction  in  the  lattice,  and  the  unit  cell  is 
of  dimensions  L\  x  L2  x  L3  =  aAi  x  (1N2  x  aN, 3, 
where  N\,  N2,  and  A3  are,  respectively,  the  num¬ 
ber  of  repeating  planes  stacked  in  the  [100]-,  [010]-, 
and  [001]-directions.  Periodic  boundary  conditions 
[74]  are  applied  along  all  faces  of  the  unit  cell  such 
that  atoms  exiting  the  unit  cell  during  the  calcula¬ 
tion  are  mapped  back  into  the  cell  on  the  opposite 
face,  thereby  preserving  the  total  mass  of  the  sys¬ 
tem.  Table  3  lists  the  details  pertaining  to  specific 


FIGURE  2.  Atomic  scale  unit  cell  for  BCC  lattice  (defect- 
free  [111]  configuration  shown) 
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TABLE  3.  Unit  cell  parameters  for  energy  minimization  computations 


Case 

# 

Defect 

Lattice 

orientation 

Atoms 

N 

Unit  cell  size 

Li  xL2  xL3  [nm] 

Vacancy 
cone,  x 

Burgers  vector 

1  b  |  [nm] 

Disl.  density 
a11  [1/nm] 

1 

vacancy 

[111] 

2016 

3.29  x  3.10  x  3.13 

5.0(10)~4 

- 

- 

2 

vacancy 

[111] 

8064 

3.29  x  6.20  x  6.27 

1.2(10)~4 

- 

- 

3 

vacancy 

[111] 

18144 

3.29  x  9.30  x  9.40 

s^ior5 

- 

- 

4 

vacancy 

[111] 

32256 

3.29  x  12.4  x  12.5 

3.1(10)”5 

- 

- 

5 

dislocation 

[111] 

8064 

3.29  x  6.20  x  6.27 

- 

0.27411 

0.0071 

6 

dislocation 

[111] 

18144 

3.29  x  9.30  x  9.40 

- 

0.27411 

0.0031 

7 

dislocation 

[111] 

32256 

3.29  x  12.4  x  12.5 

- 

0.27411 

0.0018 

8 

dislocation 

[100] 

2000 

3.17  x  3.17  x  3.17 

- 

0.31652 

0.0316 

9 

dislocation 

[100] 

8000 

3.17  x  6.33  x  6.33 

- 

0.31652 

0.0079 

10 

dislocation 

[100] 

18000 

3.17  x  9.50  x  9.50 

- 

0.31652 

0.0035 

11 

dislocation 

[100] 

32000 

3.17  x  12.7  x  12.7 

- 

0.31652 

0.0020 

unit  cells  investigated  in  this  work  and  discussed  in 
what  follows. 

First  considered  are  vacancies.  The  correspond¬ 
ing  initial  configuration  is  constructed  simply  by  re¬ 
moving  the  atom  closest  to  the  centroid  of  the  unit 
cell.  The  defect  density  in  this  case  is  defined  as  the 
volume  fraction  of  missing  atoms,  that  is,  X  =  1/N, 
where  N  is  the  total  number  of  atoms  prior  to  va¬ 
cancy  creation,  ranging  from  2016  to  32,256  among 
the  simulations  listed  in  Table  3.  Defect  energy  is 
shown  in  Fig.  3,  defined  on  a  per-atom  basis  as 


0.0020 


^  0.0004 
0.0002 

0.0000 


AEH,  x=  1 .2(1 0)*4 
CGM,  1.2(10)  4 
AEH,  x=  5.5(1 0)"5 
CGM,  5.5(1 0p 
AEH,  %=  3.1(10)"5 
CGM,  x=  3.1(1 0p 


-UHEF-O-CHEI— □— □’Q'Q—D— 


0.000  0.005  0.010  0.015  0.020  0.025 

Fu-I 


FIGURE  3.  Defect  energy  versus  applied  stretch  for  unit 
cells  with  vacancy  concentration  Energies  computed 
using  current  AEH  method  and  CGM  method  [49] 


Ed 


E-E 

N 


(59) 


where  E  is  the  total  potential  energy  of  the  system 
from  Eq.  (47)  and  E  is  the  total  potential  energy 
of  a  perfect  BCC  W  lattice  of  the  same  dimensions 
and  same  number  of  atoms  (prior  to  vacancy  for¬ 
mation),  subjected  to  the  same  deformation  bound¬ 
ary  conditions.  From  Fig.  3,  we  see  that  the  AEH 
approach  is  validated  in  the  sense  that  it  predicts 
minimum  energy  atomic  configurations  that  com¬ 
pare  favorably  with  those  obtained  using  incremen¬ 
tal  conjugate  gradient  minimization  (CGM).  Com¬ 
pared  in  Table  4  are  the  total  vacancy  energies  NEd 
at  an  applied  deformation  of  Fn  =  1.025  (2.5%  uni¬ 
axial  strain),  computed  by  AEH,  CGM,  and  the  CB 
rule.  In  these  calculations,  AEH  predicted  the  low¬ 
est  energy,  followed  by  CGM,  with  CB  predicting 
the  highest  energy.  Recall  that  all  three  methods 
commenced  from  the  same  set  of  initial  atomic  coor¬ 
dinates,  found  via  energy  minimization  from  CGM 
at  Fn  =  1.000.  For  the  CB,  the  atoms  were  dis¬ 
placed  homogeneously,  via  Fu,  from  their  initial  co¬ 
ordinates  without  any  energy  minimization  or  cor¬ 
rection,  leading  naturally  to  a  higher  defect  energy 
than  was  achieved  via  the  other  two  methods.  No¬ 
tice  from  Fig.  3  that  the  defect  energy  increases  with 
applied  stretch  Hu.  Furthermore,  from  Table  4,  the 
energy  per  vacancy  under  applied  deformation  de¬ 
creases  with  decreasing  defect  density  %.  Schultz 
[82]  reported  an  experimental  value  of  3.6  eV  for  va¬ 
cancy  formation  energy  in  pure  W  at  null  applied 
strain. 
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TABLE  4.  Vacancy  energies  at  applied  stretch  Fn  =  1.025,  computed  by  asymptotic  expansion  homogenization 
(AEH),  conjugate  gradient  minimization  (CGM),  and  the  Cauchy-Born  approximation  (CB) 


Case  # 

Defect  density  y 

NEd[ev],  Fn  =  1.025 

AEH 

CGM 

CB 

1 

5.0(10)~4 

3.6603 

3.7046 

3.7203 

2 

1.2(10)-4 

3.6492 

3.7032 

3.7189 

3 

5.5(10)-b 

3.6456 

3.7030 

3.7187 

4 

3.1(10)-5 

3.6439 

3.7029 

3.7186 

Considered  next  is  the  energy  of  W  containing 
periodic  arrays  of  screw  dislocations.  Two  types  of 
screw  dislocations  are  modeled.  In  the  first  type 
(cases  5-7  in  Table  3),  the  dislocation  tangent  line 
and  Burgers  vector  b  are  oriented  along  the  [111]- 
direction  and  pass  through  the  centroid  of  the  unit 
cell,  with  b  =  |b|  =  \fZa/2.  This  class  of  a/2  (111) 
dislocations,  most  prevalent  in  BCC  metals  such 
as  W,  is  thought  to  dominate  strain  hardening  be¬ 
havior  in  plastic  deformation  and  has  received  the 
most  attention  in  the  molecular  mechanics  litera¬ 
ture  [74,78,80,83].  In  the  second  type  (cases  8-11 
in  Table  3),  the  tangent  line  and  Burgers  vector  are 
oriented  along  the  [100]-direction  and  pass  through 
the  centroid  of  the  unit  cell,  with  b  =  |b|  =  a. 
Such  a  (100)  dislocations,  while  of  less  general  in¬ 
terest  than  the  aforementioned  a/2  (111)  type,  have 
been  observed  in  BCC  metals  [83],  but  usually  as 
part  of  hexagonal  networks,  emerging  as  a  result  of 
attractive  junctions  between  two  a/2  (111)  disloca¬ 
tions.  Nonetheless,  the  a  (100)  dislocations  are  of 
interest  here  since  they  provide  another  class  of  de¬ 
fect  whose  behavior  may  be  used  to  verify  the  accu¬ 
racy  of  the  AEH  method,  and  may  later  be  param¬ 
eterized  in  a  continuum  setting.  In  either  case,  ini¬ 
tial  atomic  positions,  prior  to  minimization  at  null 
strain  via  CGM,  are  prescribed  via  the  usual  dis¬ 
placement  field  attributed  to  a  screw  dislocation  em¬ 
bedded  in  an  infinite  isotropic  elastic  body:  u 1  = 
bd/2n,  where  0  is  an  angular  coordinate  about  the 
axis  of  the  dislocation  line.  Note  that  the  isotropic 
solution  should  be  particularly  valid  for  tungsten  as 
single  crystalline  W  is  virtually  elastically  isotropic 

[84] .  The  scalar  dislocation  density  is  defined  as  the 
defect  line  length  per  unit  reference  volume,  here 
p  =  1  /  (L2L3),  and  the  dislocation  density  tensor 

[85]  becomes,  in  this  context, 

ocAB  =  pbAE,B  (60) 


where  £,  is  the  unit  tangent  line  in  the  reference  con¬ 
figuration.  Note  that  a  denotes  the  density  of  GNDs 
in  the  sense  of  Ashby  [62].  For  the  present  set  of  sim¬ 
ulations  of  screw  dislocations,  b  ||  £,,  and  the  only 
nonvanishing  component  of  a  is  a11  =  b/  (  L->L:>,). 

The  core  structure  for  ana/2(lll)  screw  disloca¬ 
tion,  prior  to  applied  loading,  is  shown  in  Fig.  4(a). 
In  this  two-dimensional  illustration  of  the  stacking 
sequence  of  {111}  planes,  the  relative  [111]  displace¬ 
ments  of  atoms  are  denoted  by  arrows,  with  the  size 
of  each  arrow  denoting  the  magnitude  of  the  rel¬ 
ative  displacement  [74,86].  The  results  shown  for 
the  equilibrium  undeformed  structure  of  the  dislo¬ 
cation  core  are  consistent  with  the  findings  in  earlier 
literature  [74,83],  thereby  validating  the  initial  con¬ 
ditions  used  in  the  present  set  of  simulations.  The 
following  observations  are  made:  (1)  the  differential 
displacements  exhibit  threefold  symmetry  about  the 
dislocation  core;  (2)  the  largest  displacements  oc¬ 
cur  along  the  {110}  planes,  directed  radially  away 
from  the  core,  with  decreasing  magnitude  on  in¬ 
creasing  distance  from  the  core;  and  (3)  displace¬ 
ments  of  other  atoms  in  the  system  do  not  exhibit 
reflection  symmetry  about  the  {110}  planes.  Dislo¬ 
cation  energies  per  atom,  computed  via  AEH  and 
CGM,  are  compared  in  Fig.  4(b)  for  a/2  (111)  dislo¬ 
cations  and  in  Fig.  4(c)  for  a  (100)  dislocations.  The 
defect  energy  per  atom  is  defined  as  in  Eq.  (59),  that 
is,  Ed  =  (E  —  E)  /N,  where  E  is  the  total  potential 
energy  of  the  system  and  E  is  the  total  potential  en¬ 
ergy  of  a  perfect  BCC  W  lattice  of  the  same  dimen¬ 
sions  and  same  number  of  atoms,  subjected  to  the 
same  deformation  boundary  conditions.  This  quan¬ 
tity  is  computed  accurately  by  AEH,  as  is  verified 
by  close  agreement  with  the  incremental  conjugate 
gradient  (CGM)  solutions,  as  shown  in  Figs.  4(b) 
and  4(c).  Note,  however,  that  AEH  appears  to  per¬ 
form  better  (i.e.,  yields  a  lower  result  for  Ed  rela- 
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FIGURE  4.  Screw  dislocation  results,  a)  [111]  differential  displacement  map  for  a/2  (111)  screw  dislocation  [74,86], 
b)  energy  comparison  for  a/2  (111)  screw  dislocation;  c)  energy  comparison  for  a  (100)  screw  dislocation 


five  to  that  predicted  by  CGM)  for  the  a  (100)  dislo¬ 
cations  than  the  a/2  (111)  dislocations  at  large  de¬ 
fect  densities.  Furthermore,  for  a/2  (111)  disloca¬ 
tions  (Fig.  4(b)),  a  linear  increase  in  stored  defect 
energy  with  applied  deformation  Fn  is  observed. 
Flowever,  for  a  (100)  dislocations  (Fig.  4(c)),  the  de¬ 
fect  energy  remains  virtually  constant  or  decreases 
very  slightly  with  applied  deformation.  For  both 
classes  of  dislocations,  an  approximately  linear  in¬ 


crease  in  Ed  with  increasing  dislocation  density  a11 
is  evident.  Also,  comparing  Figs.  4(b)  and  4(c),  for 
similar  magnitudes  of  dislocation  density,  the  defect 
energy  for  the  a  (100)  dislocations  is  substantially 
greater,  presumably  due  to  a  larger  Burgers  vector. 
For  example,  Ed  =  0.0184  eV  /  atom  for  a  (100)  dis¬ 
locations  at  a11  =  0.0035/nm  and  at  Fu  =  1.000, 
while  Ed  =  0.0114  eV  /  atom  for  a/2  (111)  disloca¬ 
tions  at  a11  =  0.0031/nm  and  at  Fn  =  1.000. 
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5.  CONTINUUM  MODELING  OF  ENERGETICS 
OF  POINT  AND  LINE  DEFECTS 

Considered  next  are  continuum  representations  of 
energy  density  for  W  single  crystals  containing  peri¬ 
odically  distributed  vacancies  or  screw  dislocations. 
Note  that  no  attempt  is  made  to  consider  isolated 
dislocation  energies,  as  has  been  the  goal  of  previ¬ 
ous  studies  [33].  Instead,  the  effects  of  neighboring 
defects  are  purposely  included,  thereby  represent¬ 
ing  a  distribution  of  lattice  imperfections  as  would 
occur,  for  example,  in  a  plastically  deforming  sam¬ 
ple  of  material. 

In  the  context  of  continuum  defect  field  theories 
[66,87-89],  one  may  decompose  the  tangent  map 
from  an  initial,  defect-free  state  B  to  the  current  con¬ 
figuration  B  multiplicatively,  that  is, 

A  =  FK  =  FaAKAa  (61) 

where  F  is  the  compatible  deformation  gradient 
(Eq.  (1))  from  undeformed,  but  possibly  defective, 
state  Bq  to  deformed  configuration  B;  K  represents 
insertion  of  defects  into  the  macroscopically  unde¬ 
formed  lattice;  and  A  is  the  total  deformation  map¬ 
ping.  Greek  indices  denote  components  referred  to 
basis  vectors  in  defect-free  configuration  B.  The 
physics  of  Eq.  (61)  are  illustrated  conceptually  in 
Fig.  5.  Generally,  neither  A  nor  K  is  a  compatible  de¬ 
formation  mapping  as  defects  such  as  dislocations 
and  vacancies  introduce  discontinuities  in  the  lat¬ 
tice,  and  there  is  no  one-to-one  correspondence  be¬ 
tween  atoms  in  B  and  Bq.  As  such,  skew-symmetric 
spatial  gradients  of  anholonomic  mappings  A  and  K 
are  generally  nonvanishing  [54,66]. 


FIGURE  5.  Configurations  and  tangent  mappings  in 
continuum  defect  theory 


For  a  material  with  uniform  and  dilute  vacancy 
density  contraction  of  the  lattice  attributed  to  va¬ 
cancy  formation  is  of  the  isotropic  form  [54,90] 

A'l=  (l  +  xr1/36:l  (62) 

The  corresponding  strain  energy  density  per  atom, 
under  isothermal  conditions,  can  be  written  as 

*  =  *  o  +  (1  -  X)  *  e  +  X'l'v  (63) 

where  T  o  is  the  cohesive  energy  of  the  undeformed, 
defect-free  lattice;  Te  is  the  recoverable  elastic  strain 
energy;  and  Ty  is  the  formation  energy  per  va¬ 
cancy.  Note  that  T  in  Eq.  (63)  may  be  converted  to  a 
per-unit-reference-volume  basis  via  division  by  the 
atomic  volume  per  atom  (3.  To  first  order, 

vkE  =  \tAB  (. FaA  -  baA)  ( FbB  -  8bB) 

=  \<CabcdEabEcd  (64) 

where  the  elastic  constants  CAB  or  CABCD  are  eval¬ 
uated  for  a  defect-free  crystal  at  FaA  =  baA  or  Fa  b  = 
0.  Total  energies  relative  to  To,  computed  via  the 
CDF  approach  of  Eq.  (63),  are  compared  quantita¬ 
tively  with  those  obtained  from  AEH  calculations  in 
Table  5  and  Fig.  6.  For  W,  the  energetic  parameters 
entering  Eq.  (63)  are  T0  =  Ec  =  —8.90  eV  /  atom 
and  Ty  =  3.63  eV.  The  elastic  energy  density  T e 
used  in  the  CDF  approximation  is  determined  here 
by  subtracting  To  from  the  total  energy  of  the  de¬ 
formed  perfect  lattice  such  that  Eq.  (64)  is  not  used 
explicitly.  As  is  clear  from  Table  5  and  Fig.  6,  for 
simulations  1-4,  agreement  between  AEH  and  CDF 
is  within  0.1%  for  the  total  energy  T  —  To,  even  at 
the  largest  defect  density  of  x  =  5-1  (10)-4.  From 
Fig.  6,  the  total  energy  increases  with  the  applied 
deformation  in  a  quadratic  fashion,  with  initial  dif¬ 
ferences  among  the  curves  at  Fu  =  1.000  caused  by 
differences  in  initial  vacancy  concentration,  as  accu¬ 
rately  captured  by  the  third  term  on  the  right-hand 
side  of  Eq.  (63).  Defect  energy  Ed  is  listed  in  Ta¬ 
ble  6.  For  the  AEH  method,  this  is  computed  via 
Eq.  (59),  while  for  the  continuum  defect  theory,  it  is 
defined  as  summed  magnitude  of  contributions  of  % 
in  Eq.  (63)  to  elastic  and  initial  energies: 

Ed  =  X(^E  +  'fv)  (65) 

As  is  clear  from  Table  6,  defect  energies  predicted 
by  simulation  and  continuum  approximation  agree 
to  within  1%,  even  at  the  highest  vacancy  concen¬ 
tration. 
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TABLE  5.  Total  energies  described  by  multiscale  theories:  asymptotic  expansion  homogenization  (AEH,  computa¬ 
tion)  and  continuum  defect  field  (CDF,  parameter  fit) 


Case  # 

Defect 

Defect  density 

'T  —  4'o[ev/atom],  Fn  =  1.000 

T  -  4'o[ev/atom],  Fu  =  1.025 

AEH 

CDF 

AEH 

CDF 

1 

vacancy 

5.0(10)~4 

1.8021(10)-3 

1.8006(10)”3 

0.016627 

0.016604 

2 

vacancy 

1.2(10)-4 

4.5035(10)”4 

4.5015(10)”4 

0.015263 

0.015259 

3 

vacancy 

5.5(10)“5 

2.0014(10)-4 

2.0007(10)~4 

0.015012 

0.015001 

4 

vacancy 

3.1(10)~b 

1.1258(10)”4 

1.1254(10)-4 

0.014924 

0.014923 

5 

[111]  disl. 

0.0071 /nm 

0.021664 

0.022616 

0.040209 

0.038572 

6 

[Ill]  disl. 

0.0031 /nm 

0.011379 

0.010052 

0.028345 

0.025371 

7 

[Ill]  disl. 

0.0018/nm 

7.1017(10)~3 

5.6541(10)“3 

0.023323 

0.020751 

8 

[100]  disl. 

0.0316/nm 

0.114262 

0.150512 

0.131972 

0.169068 

9 

[100]  disl. 

0.0079/nm 

0.036432 

0.037628 

0.054277 

0.056184 

10 

[100]  disl. 

0.0035/nm 

0.018391 

0.016724 

0.036445 

0.035279 

11 

[100]  disl. 

0.0020/nm 

0.011235 

0.009407 

0.029438 

0.027963 

Fw 1 

FIGURE  6.  Computed  (AEH)  and  continuum  approxi¬ 
mations  (Eq.  (63))  for  total  system  energy  of  BCC  tungsten 
with  periodic  vacancy  density 


In  the  continuum  defect  field  theories,  disloca¬ 
tions  are  considered  to  be  continuously  distributed, 
and  the  strain  fields  and  displacement  discontinu¬ 
ities  of  individual  defects  are  not  represented  ex¬ 
plicitly.  Consider  a  crystal  containing  periodically 
spaced  screw  dislocations  oriented  along  the  X1- 
axis.  The  deformation  mapping  attributed  to  such 
defects  is  of  the  form 


where  the  antiplane  deformation  components  y\ 
and  y\  are  differentiable  but  generally  noninte- 
grable  functions  of  their  arguments,  and  where  G  .4 
and  g“  denote  basis  vectors  on  Bq  and  B,  respec¬ 
tively.  The  density  of  GNDs  may  then  be  computed 
from  gradients  of  K  as  follows  [50,51,55,66,91].  The 
total  Burgers  vector  B  arising  from  all  dislocations 
piercing  area  A  can  be  expressed  as  the  line  or  area 
integral 


B° 


£  K~2adXA  =  J  tABC  K~f%NcdA  (67) 

C  A 


=  /  K~}0C  aAGNcdA= 


I< 


—  lex 


NcdA 


where  C  is  a  closed  loop  encircling  A  with  unit  nor¬ 
mal  N  and  a  is  the  density  of  dislocations,  defined 
in  a  discrete  manner  in  Eq.  (60)  for  one  family  of 
dislocations  and  generalized  here  as  a  summation 
over  i  families  of  dislocation  segments  each  having 
constant  Burgers  vector  b  and  tangent  line  £,,  both 
referred  to  reference  configuration  Bq.  The  permu¬ 
tation  symbols  are  written  eABC .  From  Eq.  (67)  and 
the  identity  <9x  (KK  1 )  =  0,  the  GND  density  ten¬ 
sor  is  then 


~AC  _  KA.DBC  K-loc  _  K-\ oc.DBCKA 


(68) 


k  =  6a«ga  ®  r + y2  (*2,  a3)  g,  ®  g2 

+  Y13(X2,X3)G1®g3  (66) 


The  following  general  functional  form  of  the 
isothermal  free  energy  for  a  crystal  with  a  nonzero 
dislocation  density  is  postulated: 
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TABLE  6.  Defect  energies  described  by  multiscale  theories:  asymptotic  expansion  homogenization  (AEH,  computa¬ 
tion)  and  continuum  defect  field  (CDF,  parameter  fit) 


Case  # 

Defect 

Defect  density 

Ed[ev  /  atom 

,  Fn  =  1.000 

.Ed[ev/atom 

,  Fn  =  1.025 

AEH 

CDF 

AEH 

CDF 

1 

vacancy 

5.0(10)~4 

1.8021(10)-3 

1.8006(10)~3 

1.8156(10)~3 

1.8079(10)-a 

2 

vacancy 

1.2(10)-4 

4.5035(10)”4 

4.5015(10)“4 

4.5252(10)~4 

4.5199(10)"4 

3 

vacancy 

5.5(10)”5 

2.0014(10)-4 

2.0007(10)"4 

2.0093(10)“4 

2.0088(10)-4 

4 

vacancy 

3.1(10)~5 

1.1258(10)”4 

1.1254(10)”4 

1.1299(10)”4 

1.1300(10)-4 

5 

[111]  disk 

0.0071/nm 

0.021664 

0.022616 

0.025399 

0.023761 

6 

[111]  disk 

0.0031/nm 

0.011379 

0.010052 

0.013534 

0.010561 

7 

[111]  disk 

0.0018/nm 

0.007102 

0.005654 

0.008513 

0.005940 

8 

[100]  disk 

0.0316/nm 

0.114262 

0.150512 

0.113417 

0.150512 

9 

[100]  disk 

0.0079/nm 

0.036432 

0.037628 

0.035721 

0.037628 

10 

[100]  disk 

0.0035/nm 

0.018391 

0.016724 

0.017889 

0.016724 

11 

[100]  disk 

0.0020/nm 

0.011235 

0.009407 

0.010882 

0.009407 

^  =  ^o  +  'fE  +  ^{l2ocABGACGBDOcCDy  (69) 

where  G  is  a  (dimensionless)  metric  tensor  defined 
on  B0  used  to  compute  the  inner  product  of  the 
contravariant  tensor  a,  p  is  an  elastic  shear  mod¬ 
ulus  (constant),  and  l  (units  of  length)  and  m  (di¬ 
mensionless)  are  scalars.  For  W  of  interest  here, 
li  =  159  GPa  =  16.07  eV  /  atom.  Typically,  the  met¬ 
ric  in  Eq.  (69)  is  selected  for  simplicity  as  Gab  = 
8ab  [51-53];  however,  there  is  some  evidence  that 
other  choices  may  be  more  appropriate.  Gibeling 
and  Nix  [92]  discussed,  from  the  standpoint  of  dis¬ 
crete  dislocation  modeling,  how  the  strain  energy 
sustained  by  dislocations  may  be  amplified  by  ex¬ 
ternally  applied  deformations.  As  such,  following 
[66,93],  the  covariant  elastic  deformation  measure 
Gab  =  Gab  =  FaAgabFbB  may  be  used  in  Eq.  (69)  to 
reflect  amplification  of  internal  defect  energy  com¬ 
mensurate  with  applied  deformation.  The  choice 
of  m  in  Eq.  (69)  also  warrants  careful  consideration. 
Most  often  in  recent  literature,  m  =  1  is  used  to  re¬ 
flect  a  quadratic  dependence  of  the  free  energy  on 
the  GND  tensor  [51-53].  This  form  is  often  invoked, 
in  conjunction  with  thermodynamic  arguments,  to 
provide  a  back  stress  dependent  on  the  density  ten¬ 
sor  of  dislocations  [52,64]  and  /  or  its  spatial  gradient 
[51,53].  Clayton  [63]  assumed  a  linear  dependence 
(m  =  1/2)  of  stored  energy  on  dislocation  density  in 
a  continuum  crystal  plasticity  model  of  single  crys¬ 
talline  W.  A  linear  dependence  of  free  energy  on 
SSD  density  is  also  commonly  prescribed  [52].  Fi¬ 


nally,  in  continuum  dislocation-based  plasticity  the¬ 
ories,  the  length  parameter  l  is  often  calibrated  in  a 
phenomenological  manner  to  reflect  the  degree  of 
additional  stiffness  or  strain  hardening  imparted  by 
the  dislocation  density,  for  example,  in  torsion  of  a 
thin  wire  [94]  or  in  nanoindentation  [95,96]. 

In  the  present  work,  parameters  l  and  m  enter¬ 
ing  Eq.  (69)  are  computed  directly  via  a  best  fit  to 
results  from  AEFI  calculations,  as  opposed  to  a  cali¬ 
bration  to  macroscopic  data.  For  both  a/2  (111)  and 
a  (100)  dislocations,  m  =  1/2  provides  a  superior  fit 
to  the  computational  results  than  does  m  =  1.  Since 
the  defect  energy  Ed  in  Fig.  4(b)  for  a/2  (111)  dislo¬ 
cations  increases  with  applied  deformation.  Gab  = 
Cab  =  FaAgabFbB  is  used  for  dislocations  of  that  ori¬ 
entation.  Note  that  for  uniaxial  strain  of  the  form 
Fn  =  1  +  £,  with  screw  dislocations  oriented  paral¬ 
lel  to  the  A1 -direction,  the  defect  energy  in  Eq.  (69) 
degenerates  to 

Ed  =  p  (i2ocabGAcGbdoccd)  1/2 

=  [ila11  (l  +  21  +  e2)  (70) 

On  the  other  hand,  since  the  energy  Ed  in  Fig.  4(c) 
of  a  (100)  dislocations  remains  relatively  constant 
with  applied  loading.  Gab  =  &ab  is  used  for  rep¬ 
resenting  the  energy  of  a  (100)  dislocations,  mean¬ 
ing  that  Eq.  (69)  applies  in  that  case  with  l  =  0. 
For  a/2  (111)  dislocations,  l  =  0.200  nm  =  0.736, 
and  for  a  (100)  dislocations,  l  =  0.297  nm  =  0.946. 
These  values  are  much  smaller  in  magnitude  than 


Volume  5,  Number  3&4,  2007 


220 

Begell  House  Inc.,  http://begellhouse.com  Downloaded  2008-1-29  from  IP  128.63.66.91  by  Dr.  John  Clayton  (jdclayt) 


CHUNG  AND  CLAYTON 


some  of  those  calibrated  from  macroscopic  data, 
with  l  reported  in  the  latter  on  the  order  of  several 
micrometers  [94,95],  though  values  in  the  nanome¬ 
ter  range  have  been  suggested  for  l  from  hardness 
measurements  inferred  from  indentation  data  [96]. 
Comparisons  of  the  energy  4'  —  't'0  described  by 
the  CDF  theory  of  Eq.  (69)  with  numerical  predic¬ 
tions  from  AEH  over  the  range  of  uniaxial  deforma¬ 
tion  1.000  <  F\  i  <  1.025  are  shown  in  Fig.  7  for 
a/2  (111)  dislocations  and  in  Fig.  8  for  a  (100)  dislo¬ 
cations.  For  ease  of  quantitative  comparison,  values 
at  Fn  =  1.000  and  Fn  =  1.025  are  tabulated  for 
the  relative  total  energy  4'  —  4' o  in  Table  5  and  de¬ 
fect  energy  Ed  in  Table  6.  The  agreement  between 
AEH  and  CDF  predictions  is  generally  modest  and 
is  closest  at  relatively  large  defect  densities,  that  is, 
for  a11  >  0.005  /  nm. 

The  investigation  conducted  here  is  limited  in 
the  sense  that  only  a  few  classes  of  screw  dislo¬ 
cations  are  examined,  for  only  one  material,  and 
these  defects  are  arranged  periodically  in  the  lat¬ 
tice.  Furthermore,  all  simulations  are  isothermal  at 
null  temperature  (akin  to  molecular  statics).  How¬ 
ever,  the  work  represents  an  initial  step  toward  com¬ 
puting  energies  used  in  continuum  defect  theories 
from  physics-based,  multiscale-atomistic  computa¬ 
tions,  as  opposed  to  phenomenological  curve  fit- 


Cu-1 

FIGURE  7.  Computed  (AEH)  and  continuum  approx¬ 
imations  (Eq.  (69))  of  energy  for  BCC  tungsten  with  non¬ 
zero  dislocation  density  component  a11  arising  from  pe¬ 
riodic  a/2  (111)  screw  dislocations 
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FIGURE  8.  Computed  (AEH)  and  continuum  approx¬ 
imations  (Eq.  (69))  of  energy  for  BCC  tungsten  with  non¬ 
zero  dislocation  density  component  a11  arising  from  pe¬ 
riodic  a  (100)  screw  dislocations 


ting  to  macroscopic  stress-strain  or  microindenta¬ 
tion  data,  for  example.  The  dislocation  densities 
considered  here,  with  p  on  the  order  of  1016/m2, 
are  quite  large  relative  to  densities  found  in  homo¬ 
geneously  deforming  single  crystals.  Argon  and 
Maloof  [97]  reported  values  of  p  on  the  order  of 
1012/m2  for  pure  tungsten  single  crystals  deformed 
up  to  5%  axial  tensile  strain.  Densities  examined 
presently  correspond  to  a  dislocation  spacing  on 
the  order  of  10  nm.  In  highly  strained  regions  of 
crystal,  such  as  in  the  vicinity  of  grain  or  subgrain 
boundaries  formed  during  severe  plastic  deforma¬ 
tion  (SPD),  such  tight  spacing  may  not  be  unrea¬ 
sonable.  For  example,  if  one  considers  a  bound¬ 
ary  comprising  a  sequence  of  dislocations  of  spac¬ 
ing  h  and  Burgers  vector  b,  the  misorientation  at 
the  boundary  can  be  computed  as  b  /h  [98],  on  the 
order  of  4°  for  10-nm-spaced  dislocations  in  tung¬ 
sten.  For  subgrain  boundaries  produced  in  tungsten 
crystals  deformed  through  SPD  processes  [99,100], 
misorientations  of  such  magnitude  have  been  doc¬ 
umented.  Furthermore,  strain  gradient-based  con¬ 
tinuum  defect  theories  are  designed  to  address  such 
phenomena  in  the  context  of  crystal  plasticity  [51- 
55],  so  results  presented  here  may  be  applied  to  mo¬ 
tivate  continuum  energy  dependencies,  at  least  in 
a  qualitative  sense,  on  defect  densities  that  serve 
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an  important  role  in  such  theories.  The  thermody¬ 
namics  of  stored  energy  of  cold  working  may  also 
influence  shear  localization  processes  in  ultra  fine 
grained  tungsten  [100,101],  a  material  that  can  ex¬ 
hibit  dislocation  densities  of  the  magnitudes  stud¬ 
ied  here. 

6.  CONCLUSIONS 

A  multiscale  method  based  on  the  theory  of  asymp¬ 
totic  expansion  homogenization  is  developed  and 
implemented.  The  technique  enables  computation 
of  global  or  coarse-grained  mechanical  properties 
such  as  effective  elastic  stiffness  and  net  stress  of 
a  deforming  unit  cell  consisting  of  periodically  ar¬ 
ranged,  discrete  atoms  at  the  fine  scale.  The  for¬ 
mulation  directly  accounts  for  the  effects  of  defects 
in  the  lattice.  From  a  computational  standpoint, 
the  method  is  an  efficient  alternative  to  incremen¬ 
tal  conjugate  gradient  schemes  in  terms  of  predic¬ 
tion  of  minimum  energy  configurations  of  atomic 
degrees  of  freedom  in  statically  deforming  lattices 
containing  defects.  The  direction  of  a  perturba¬ 
tive  displacement — which  serves  as  a  correction  to 
the  CB  rule  which  is  accurate  only  for  homoge¬ 
neous  deformations — is  computed  via  application 
of  the  homogenization  scheme  in  an  equation  that 
relates  the  free  energy  of  the  system,  the  macro¬ 
scopic  deformation  gradient,  and  local  atomic  de¬ 
grees  of  freedom.  The  key  computational  advan¬ 
tage  of  this  approach,  relative  to  traditional  molecu¬ 
lar  statics  or  molecular  dynamics,  is  the  use  of  a  di¬ 
rected  line  search  algorithm,  which  reduces  the  en¬ 
ergy  minimization  process  from  minimization  over 
traditional  3  x  N  solution  space  to  a  line  search  in 
one  dimension  (with  this  direction  determined  as 
mentioned  above  via  solution  of  the  fine-scale  equi¬ 
librium  equation  of  the  homogenization  scheme). 
The  current  multiscale  method  also  does  not  present 
any  difficulties  associated  with  matching  boundary 
conditions  across  length  scales  since  the  displace¬ 
ment  gradient,  rather  than  displacements  them¬ 
selves,  is  the  kinematic  quantity  exchanged  between 
scales.  In  validation  simulations  with  predefined 
defect  structures,  close  agreement  is  found  with  en¬ 
ergies  predicted  by  the  present  method  and  conju¬ 
gate  gradient-based  molecular  mechanics,  and  cat¬ 
egorically  lower  system  energies  are  predicted  than 
those  obtained  from  the  CB  rule. 

Specifically  studied  here  are  the  nonlinear  elas¬ 
tic  responses  of  BCC  tungsten  single  crystals  con¬ 


taining  periodically  distributed  vacancies  and  screw 
dislocations  of  two  different  classes.  It  is  found  that 
energies  associated  with  vacancies  and  a/2  (111) 
screw  dislocations  tend  to  increase  with  applied 
uniaxial  stretching,  while  energies  of  a  (100)  screw 
dislocations  tend  to  remain  constant  with  stretch. 
These  computed  energies  are  used  to  motivate  con¬ 
tinuum  energy  functions  for  defective  crystals  de¬ 
scribed  in  terms  of  the  vacancy  density,  the  disloca¬ 
tion  density  tensor,  and/or  the  applied  deformation 
gradient.  For  crystals  with  vacancies,  a  continuum 
description  of  defect  energy  increasing  linearly  with 
vacancy  density  is  found  to  be  extremely  accurate, 
relative  to  the  computational  results,  over  the  range 
of  defect  densities  considered.  For  W  with  screw 
dislocations,  a  defect  energy  linearly  dependent  on 
dislocation  density  provides  better  agreement  with 
numerical  results  than  does  a  quadratic  dependency 
of  this  energy  on  the  dislocation  density  tensor  often 
encountered  in  dislocation-based  gradient  plasticity 
theory. 

The  method  and  results  presented  here  offer  sev¬ 
eral  avenues  for  improvement  of  existing  crystal 
elasticity  and  plasticity  models.  The  present  fo¬ 
cus  has  been  on  the  nonlinear  elastic  response  as 
well  as  construction  of  elastic  strain-  and  defect- 
energy  functions  entering  continuum  gradient  plas¬ 
ticity  theories.  The  results  could  be  used  in  stand¬ 
alone  crystal  plasticity  theory,  either  classical  local 
theory  [102]  or  nonlocal  gradient  theory  [55],  in  ei¬ 
ther  case  with  the  stored  energy  and  elastic  moduli 
allotted  an  explicit  dependence  on  initial  defect  con¬ 
tent.  More  ambitiously,  the  present  homogeniza¬ 
tion  technique  could  be  incorporated  in  a  discrete- 
continuum  multiscale  context,  with  elastic  modu¬ 
lus  and  stored  energy  depending  on  the  fine-scale 
configuration  of  defects  resolved  concurrently  in  the 
simulation;  such  an  idea  was  pursued  in  [45]  for  sin¬ 
gle  slip  conditions.  Finite  temperature  effects  would 
be  needed  for  realistic  resolution  of  dislocation  dy¬ 
namics.  Incorporation  of  such  effects  would  involve 
some  reformulation  of  the  governing  equations  to 
account  for  dynamic  as  opposed  to  static  conditions, 
beginning  with  Eq.  (4),  or  consideration  of  alterna¬ 
tive  computational  strategies  for  irreversible  ther¬ 
modynamics,  as  discussed  below.  In  principle,  ther¬ 
mally  activated  dislocation  line  motion  could  then 
be  transferred  upward  in  scale  to  describe  contin¬ 
uum  plastic  deformation  [50]. 
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Finally,  it  is  appropriate  to  note  the  scope  of 
potential  applications  and  some  of  the  limitations 
of  the  method  proposed  in  this  article.  Consider 
specifically  the  developments  leading  to  Eq.  (32) 
and  several  of  the  general  underlying  assumptions 
of  the  homogenization  technique.  The  use  of  a  trun¬ 
cated  expansion  of  the  strain  energy  presumes  that 
the  system  stays  in  the  positive  definite  region  of  the 
energy  landscape  and  that  the  deformations  con¬ 
sidered  are  primarily  of  the  thermodynamically  re¬ 
versible  kind.  While  it  is  possible  to  incorporate 
irreversible  mechanisms,  as  has  been  attempted  in 
[45]  by  prescribing  specific  kinetic  pathways  for  dis¬ 
location  glide,  the  development  of  general  models 
accounting  for  irreversible  continuum-scale  mecha¬ 
nisms  such  as,  say,  plasticity,  based  on  atomistic  in¬ 
teractions,  has  yet  to  be  fully  realized.  Such  mod¬ 
els  would  be  needed  to  account  for  dissipation  and 
the  nonequilibrium  effects  associated  with  energy 
transfer  in  and  around  the  simulation  domain.  In 
such  cases,  modeling  approaches  using  a  gener¬ 
alized  homogenization  method  to  account  for  the 
time  variable  may  be  needed  [46]. 

The  limitations  imposed  on  possible  applications 
by  the  assumption  of  periodicity  merit  reiteration. 
In  this  work,  the  primary  results  stem  from  the  ac¬ 
counting  of  the  physically  inhomogeneous  defor¬ 
mation  around  crystal  defects,  with  a  demonstrated 
connection  between  atomic  and  continuum  model¬ 
ing  methods.  Insofar  as  the  distribution  of  imper¬ 
fections  such  as  dislocations  can  be  described  as  be¬ 
ing  periodic,  the  proposed  homogenization  method 
is  applicable  [33],  with  suitable  accuracy,  as  sug¬ 
gested  by  the  results  reported  here.  In  the  event 
that  the  system  departs  from  periodicity,  a  common 
occurrence  in  the  irreversible  regime,  where  dislo¬ 
cations  move  and  can  ultimately  pile  up  along  in¬ 
terfaces,  adaptive  modeling  approaches  [103-106] 
may  be  most  suitable.  In  such  instances,  hetero¬ 
geneous  solution  domains  involving  some  subdo¬ 
mains  modeled  through  homogenization,  as  de¬ 
scribed  here,  and  others  modeled  with  quasicon¬ 
tinuum  theory  [26-28]  or  direct  molecular  statics, 
for  example,  would  be  appropriate.  Investigations 
of  the  like  suggest  interesting  directions  for  future 
work.  With  the  present  homogenization  technique, 
it  should  also  be  noted  that  periodicity  of  the  lat¬ 
tice  enables  one  to  consider  a  reduced  number  of 
atoms,  relative  to  a  more  disordered  material,  at  the 
fine  scale.  For  example,  modeling  of  a  single  crys¬ 


tal  of  uniform  lattice  orientation  generally  requires 
far  fewer  atoms  than  would  realistically  modeling 
a  polycrystal.  The  latter,  though  conceptually  fea¬ 
sible,  would  require  computationally  costly  atom¬ 
istic  resolution  of  grain  boundaries  and  intergranu¬ 
lar  misorientations. 
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